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ABSTRACT 

Many galaxy clusters pose a "cooling-flow problem" , where the observed X-ray emis- 
sion from their cores is not accompanied by enough cold gas or star formation. A 
continuous energy source is required to balance the cooling rate over the whole core 
volume. We address the feasibility of a gravitational heating mechanism, utilizing the 
gravitational energy released by the gas that streams into the potential well of the 
cluster dark-matter halo. We focus here on a specific form of gravitational heating in 
which the energy is transferred to the medium thorough the drag exerted on inflowing 
gas clumps. Using sphcri-symmctric hydro simulations with a subgrid representation 
of these clumps, we confirm our earlier estimates that in haloes ^ 10 13 M Q the grav- 
itational heating is more efficient than the cooling everywhere. The worry was that 
this could overheat the core and generate an instability that might push it away from 
equilibrium. However, we find that the overheating does not change the global halo 
properties, and that convection can stabilize the cluster by carrying energy away from 
the overheated core. In a typical rich cluster of 10 14_15 M Q , with ~ 5% of the accreted 
baryons in gas clumps of ~ 10 8 M Q , we derive upper and lower limits for the tem- 
perature and entropy profiles and show that they are consistent with those observed 
in cool-core clusters. We predict the density and mass of cold gas and the level of 
turbulence driven by the clump accretion. We conclude that gravitational heating is 
a feasible mechanism for preventing cooling flows in clusters. 

Key words: galaxies: clusters: general — galaxies: haloes — (galaxies:) cooling flows 
— galaxies: formation — hydrodynamics — X-rays: galaxies: clusters 
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1 INTRODUCTION 

Galaxy clusters can be divided into two distinct popula- 
tions according to the X-ray luminosity of their central 
cores ([Sanderson et al.ll2006[ ). Cool-core (CC) clusters 
are centrally concentrated, highly luminous in X-ray, 
and have central cooling times of 0.1 — 1 Gyr. Non- 
cool-core (NCC) clusters have lower densities and lu- 
minosities near the centre, and their central coolin g 
times are typically a few Gyrs (jDonahue et al.|[2006t ). 
The population of CCs has little internal variability, 
and they all exhibit a typical temperature profile with 
a decline by a factor of 2-3 in the innermost few 10 kpc 
(jLeccardi fe Molendill2008h . Based on the short cooling 
time in CCs, one expects to observe gas in intermittent 
temperatures, a high star-formation rate in the bright- 
est central galaxy (BCG) (> 1OOM / yr), and a large 



stellar mass in the BCG (1O 12 ~ 13 M ), none of which is 
observed. These are thre e manifestati ons of the cooling- 
flow problem in clusters (|Fabianll994l) . Since the cooling 
time is inferred directly from observations based on the 
luminosity and temperature, the discrepancy between 
the expected cooling rates and the gas that actually 
cools indicates that some heating mechanism is keeping 
the gas hot in a stable configuration. 

Several mechanisms have been proposed as poten- 
tial solutions to the cooling-flow problem. Most popular 
is AGN feedback, where energy or momentum are pro- 
vided by an a ctive galactic nucleus ei ther in an intense 
quasar mode (ICiotti fc Ostrikerll2007ri. in a slower radio 
mode (jBes t 2007; Catt aneo et al.ll2009l for a review), or 
via cosmic rays (|Guo fc Ohll2008[) . The AGNs clearly re- 
lease sufficient power to balance the cooling in the cores, 
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and the obser ved big radio and X -ray bubbles in some 
cluster cores (jBirzan et al.1120041) is likely evidence for 
AGN feedback. However, the coupling of the AGN en- 
ergy to the gas in the whole core volume is not easily 
understood, and the requir ement of continuous he ating 
is not a trivial constraint (|De Young et al.ll2008f ). An- 
other possibility is that the cooling instability might 
be locally suppressed with the addition of non-thermal 
pressure so urces such as cosmic rays or turbulent mag- 
netic fields (jSharma et al.ll2010l ). 

Here we utilize the fact that the gravita- 
tional power released as fresh baryons stream into 
the potential well created by a massive halo is 
more than enough t o ba la nce the cooling rate at 
the centre (iFabianl 120031: iDekel k Birnboiml 120081: 
IWang k Abell 120081) . Among the mechanisms through 
which this energy is transferred to t he inner halo 
one c ould consider dynamical friction (lEl-Zant et al.l 
2004 iFaltenbacher k Mathewsll2007t iNaab et al 1120071; 



Khochfar k Ostrikerl 120081: iJohansson et al.l I2009D . 
therm al conduction iKim k Naravanl ( 20031 ). or turbu- 
lence (|Sato k Nagatakil I2004T ) . In particular, the tur- 
bulence induced by accreting sub-structures may dis- 
turb the magnetic field in a w ay that could make 
the conduction more efficient (jParrish et al.l 120101: 
iRuszkowski k Ohll2010D . 

The gravitational power that is released during the 
streaming of baryons into clusters of M ^ 10 13 M Q is 
indeed sufficient for balancing the expected radiative 
losses (jDekel k Birnboiml 120081 hereafter DBO80- The 
challenge is to deposit this energy at the inner cluster 
core of ~ 30 — 100 kpc, where most of the cooling takes 
place. The energy should be evenly distributed over the 
whole core volume and continuously over several Gi- 
gayears. It should be performed in a wa y that is consis- 
tent with the observe d entropy profile (|Donahue et al.l 
120061 ) and metallicity (jRebusco et al.ll2006| ) profile. Here 
we investigate a model in which some of the accreted 
gas is in dense and cold (~ 10 4 K) gas clumps. These 
clumps do not stop at the virial shock but rather pen- 
etrate to the inner parts of the halo. The continuous 
accretion of many clumps through the halo can dis- 
tribute the energy smoothly in a large volume, as re- 
quired, unlike AGN feedback that is episodic and orig- 
inates from a very small region near the black hole. 
The main physical mechanism that couples the clumps 
and the halo gas is hydrodynamic drag, subsonic or 
supersonic, which decelerates the clumps and causes 
the deposition of their kinetic and potential energy in 
the ambient medium. Hydrodynamic drag is more effi- 
cient when the clumps are smaller and moving faster, 
as opposed to dynamical friction that becomes more ef- 
ficien t for more mas sive clumps and for transonic veloc- 
ities (IOstrikerlll999h . We pointed out in DB08 that the 
gravitational heating is likely to be associated with the 
streams tha t build the cluster along the filaments of th e 
cosmic web (jDekel k Birnboiml2006HDekel et al.ll2009D . 
These streams could transfer their energy to the cluster 



1 We note that the overcooling problem is less pron ounced for 
haloes between 10 12 - 10 13 M Q jBirnboim et alj|2007l) . 



core via the generation of turbulence and other mech- 
anisms that are not neces sarily associated with cold 
clumps l|Zinger et al.l 1201 lb . Still, the heating through 
cold clumps is a concrete example that allows a simple 
study of the dynamical response with general implica- 
tions that are not limited to this particular coupling 
mechanism. 

In DB08 we showed that accreted clumps could 
balance the cooling in haloes more massive than 6 x 
10 12 M Q , provided that the clumps contain a non- 
negligible fraction of the infalling gas 10% for 
a 10 15 M Q halo), and that the clump masses are in 
the range 10 4 — 10 8 M Q . The clumps heat the intra- 
cluster medium (ICM) via drag u ntil they disintegrate 
by hydrodynamical in stabilities ([Murray k Linl 12004 
iMaller k Bullockll2004) . The permitted mass range for 
the clumps took into account additional criteria for 
clump survivability, such as Bonnor-Ebcrt instability, 
conduction and evaporation. When the clumps are suf- 
ficiently massive, they survive long enough to penetrate 
through the outer halo and reach the centre, while by 
being not too massive, their interaction with the ICM 
is sufficiently strong for most of the clump energy to be 
deposited in the core within a Hubble time. The basic 
features studied in DB08 assuming a static system will 
be implemented below in a dynamical evolving config- 
uration. 

The presnt study is motivated by a potential over- 
heating problem. The heating rate by drag (and by 
other mechanisms such as dynamical friction, cosmic 
rays from AGN) is proportional to the density of the 
hot ambient gas, ehoat oc phot- On the other hand, the 
Brcmsstrahlung cooling rate per unit volume scales like 
ecooi oc p h ^ t (assuming isobaric gas). This gen erates 
an instability through a posit ive feedback loop (jFieldl 
119651: IConrov k Ostrikerl l2008t) . A small negative den- 
sity perturbation, e.g., produced by gas expansion due 
to overheating, would make the ratio of heating to cool- 
ing rate increase as e coo i/ehcat oc p° D 5 t , leading to more 
overheating, enhanced pressure, and a runaway expan- 
sion. An analysis of the consequences of this instability, 
and an investigation of possible mechanisms that could 
keep the cluster in equilibrium, necessitate a dynamical 
treatment. 

This unstable overheating, as reproduced in spheri- 
symmetric simulations below, results in expanding 
shells that heat to temperatures as high as lO 9 !^, 
clearly at odds with observations. In the real world, 
this overheating must be regulated by processes that 
smooth temperature or entropy gradients by heat 
transfer through conduction or convection. Conduc- 
tion is suppre ssed in the presence of magnetic fields 
(jFabianl Il994 and reference within), though it might 
be boosted a lmost to its ma ximum possible value 
(ISpitzerl I1962D by turbul ence dNaravan k Medyedeyl 
200 It iBalbus k Reynolds! l2008t iParrish et al.l I201C 



Ruszkowski k Ohl 120101) . which might be a natural 



product of clump heating. Convection is a promising 
mechanism for smoothing the local instabilities. In the 
simple case of an ideal gas with a uniform chemical com- 
position in a spherical potential well, convection occurs 
in regions where the entropy is declining with radius. 
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However, the strength of this convection is uncertain, 
as it depends on the gradients in gas properties and 
on the magnetic fields. Weak magnetic fields make the 
gas more susceptible to convection, with the entropy- 
gra dient criterion replaced by the tem perature gradi- 
ent (|Balbusll2000l 120011 : lQuataerdl2008l ) , but for certain 
types of perturbations, the convection strength might 
be drasti cally suppressed by effects related to magnetic 
tension (|Parrish et al.l 120081 I2009D . Regardless of the 
actual energy transport mechanism, one expects nature 
not to permit shells of ~ 10 9 K in close contact with 
shells of ~ 10 7 K, and to act to smooth such a disconti- 
nuity. Motivated by the turbulence that is expected to 
be generated by the clumps, we focus below on convec- 
tion as the mechanism that smooths steep gradients. 

In this paper, we mimic this smo othing proces s by a 
ID mixing-length convection model (Spiegel 1963), with 
a the mixing-length coefficient L the single free param- 
eter. We find that the results are almost independent of 
the value of this parameter as long as noticeable convec- 
tion occurs. This allows us to further simplify the model 
by assuming that the convection is maximal, namely the 
energy transfer rate is limited by the requirement that 
hot bubbles accelerate until they become supersonic, at 
which point they disintegrate. This leaves us with no 
free parameters in our convection model. 

There are some additional benefits from a dynamic 
treatment of the clump heating process. The analysis 
of DB08 considered simple Monte-Carlo clump trajec- 
tories within an otherwise static halo in hydrostatic 
equilibrium. This assumption of a static halo could be 
valid for one Gigayear but the cluster may evolve con- 
siderably over a Hubble time, due, for example, to the 
gradual increase of virial temperature and the growth 
of the BCG. Furthermore, if clump heating is taking 
place, cold clumps continually get destroyed near the 
halo centre, dumping cold gas near the BCG. This di- 
lution of the hot gas with cold gas is not a problem for 
the heating-cooling balance because the clumps bring in 
several times the energy needed for heating themselves 
to the cluster virial temperature, but a proper account 
of this clump deposition requires a dynamical analysis 
of an evolving cluster. 

In [J3]wc describe the implementation of the clump 
model and of the convection model in the ID hydrody- 
namic code. In iJ3]we show results of hydrodynamic clus- 
ter simulations with convection and clump heating that 
match the observed temperature and entropy profiles of 
clusters and the cooling rates in clusters without a need 
for any additional feedback. In S|4]we address possible 
direct and indirect observations of the cold clumps. In 
|JS]we summarize and conclude. 



2 METHODS 

2.1 Implementing Clumps in ID Hydrodynamic 
Simulations 

According to our estimates in DB08, heating by clumps 
requires clump masses in the range 10 4 — 10 8 A/q. In 
order to properly resolve drag forces and clump disinte- 
gration via hydrodynamical instabilities, these clumps 



should be resolved by at least 1000 cells or SPH parti- 
cles (i.e., 10 cells across each dimension). The implied 
required dynamical range in a cluster of 10 14 — 1O 15 M 
is impractical, so simulation of such clumps requires a 
sub-grid model. We develop such a model below, and 
describe its implementation in a ID spherical hydrody- 
namical code. 

The clumps are made of cold and partly ionized 
gas at ~ 10 4 K in pressure equilibrium with their sur- 
rounding hot halo. For a rich cluster of galaxies, with a 
virial temperature ~ 10 7 K, the overdensity within the 
clumps is about 10 3 . The clumps couple to the hot gas 
by a drag force 

fdrag = ~ Cd A p hot l^el , (1) 

acting opposite to the relative direction of motion. Here 
is the drag coefficient (~ 1 for a spherical gas clump), 
A is the cross-section surface area of the clump (irR^), 
Phot is the density of the hot component, and v rc \ is 
the relative velocity between the clump and the halo 
gas. Eq. (fTJ) holds for subsonic and supersonic motions, 
though the value of Cd may vary, especially in the trans- 
sonic regime where it could be a few. The deceleration 
(fdrag /rn c i) is proportional to the ratio of clump surface 

— 1/3 

area to volume, oc m cl . This dictates lower and upper 
limits to the relevant clump masses. Clumps that are 
too small cannot penetrate through the outer halo into 
the core, and clumps that are to large cannot deposit a 
significant fraction of their energy in the inner halo on 

a time scale shorter than the Hubble time. 

Single clump simulations (jMurrav fc Lin! 120041) in- 
dicate that most of the energy dissipated in this process 
goes into the hot ambient gas. The survivability of these 
clumps is an open question as they could be destroyed 
by many diffe rent effects. This has be en discussed in 
DB08 and in iMaller fc Bullockl (|2004l) . In a nutshell, 
if the clump is more massive than ~ 10 s Mq it would 
exceed the Bonnor-Ebert critical mass (the equivalent 
of the Jeans mass for pressure-confined spheres) and 
it would collapse under its self-gravity and turn into 
stars. If the clump is less massive than ~ 10 4 M Q , evap- 
oration and conduction are expected to disintegrate the 
clump. Hydrodynamic instabilities, particularly Kelvin- 
Hclmholtz instability, tend to break the clump, typically 
after it has r epelled the equivalen t of its own mass in 
ambient gas ((Murray fc Lin! [20041 ) . The clump masses 
then cascade down toward the lower limit for clump 
mass. 



2.1.1 The Hydra Code 

Subgrid recipes of the effects mentioned above were 
incorporated into the ID spherical code Hydra 
(jBirnboim fc D ckel 200j|). The code is finite-difference 
Lagrangian with von-Ncumann first and second order 
artificial viscosity. Dark matter is described as zero- 
width shells that propagate through the gaseous shells, 
and interact with them gravitationally. The coupled 
density fields of gas and dark matter shells are prop- 
agated using a 4th order Runge-Kutta method. Time 
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steps are defined by the minimum of the Courant con- 
ditions and the allowed deviation of the forth order 
scheme from fully explicit (1st order) time step. When 
this difference exceeds some preset epsilon, the previ- 
ous values arc restored, and a new step with decreased 
time step is performed. A comprehensive description of 
the hydrodynamic and dynamic equations for the gas 
and the dark matter shells, the numerical scheme used 
and convergence tests and a comparison to analytic test 
pr oblems can be found in se ction 3 and the appendixes 
of iBirnboim fc Dekell (l2003h . 

The baryons and the dark matter shells arc assigned 
a preset angular momentum that is added as an im- 
pulse when the shell is at its turn-around. This pre- 
vents a numerical and physical singularity at the cen- 
tre. The angular momentum of the dark matter is set 
so that the rotational kinetic energy is 18% of the ra- 
dial kinetic energy at the virial radius. The results are 
insensitive to this choice. The baryonic angular momen- 
tum is set to produce a spherical, angular-momentum 
supported "disc" of radius ~ lOkpc for a cluster halo 
of 10 15 M Q , to mimic a BCG. Since centrifugal forces 



at the centre. Clump shells typically oscillate around 
the halo's centre before the processes described below 
destroy them. 

Each shell is assumed to contain ?i c i clumps with 
mass 



M 



shell 



"cl 



(2) 



Msheii and m c are the total shell mass and the mass 
of each gaseous clump respectively. The shells interact 
with the baryons by decelerating according to eq. ([1]). 
The drag force equation and energy equation of an in- 
teraction between some clump i and a parcel of gas are: 



and 



E' 



ft 



-F, 



gas • 



0, 



(•3) 



(4) 



respectively with f l cX and -Fg as the forces on the clump 
and gas parcels arising from the clump-gas interactions, 



scale like r" 3 , the angular momentum of the baryons and v l\ and u «« the velocities of clump i and the gas re- 



is negligible at a distance comparable to a few disc 
radii above the disc. We note that a spherical code is 
not the right tool for studying disc formation, and this 
setup is essentially an inner boundary condition for the 
halo simulation. In addition, the code imposes a cen- 
tral smoothing length on the gravitational acceleration, 
a g = G M/(r + s) 2 , typically with s ~ 50 pc. Radiative 
cooling is calculated by a metallicity d ependent cooling- 
function ([Sutherland fe Dopital fl993) with a constant 
preset metallicity. 

The initial conditions, in terms of shell masses, radii 
and peculiar velocities, were set at z = 100 such that the 
future accretion rate on to the grow i ng halo will follow 
a des ired accretion rate (|Dekellll98lt IBirnboim fc Dekel 
l2003l appendix C). Specifically, the initial perturbation 
used here yields an accretion history that traces that 
of an average main progenitor accordi n g to the EPS 
appro xi mationlPress fc Schcchtcr] (|1974[ ); lLacev fc Cold 
(|1993t) ; iNeistein et all (120061). The proce dure is de- 
scribed in detail in IBirnboim et ahl (|2007|) . The code 
ha s been compared su ccessfully to analytic predictions 
of iBertschingerl (|1985l ) and to a Von-Neumann-Scdov- 
Taylor problem. The Courant conditions and epsilons 
are set so that the global energy conservation over a 
Hubble time is always better than 1% and the spa- 
tial convergence was tested for each set of simulations. 
Timesteps arc consequently ~ 10 -5 Gyr throughout 
most of the simulation. 



2.1.2 Drag forces in ID 

The subgrid model for clumps is similar in its approach 
to "sticky particle" techniques in the sense that it cal- 
culates subgrid interactions on otherwise N-body par- 
ticles. Here we define "clump-shells", which, like dark 
matter shells, arc able to penetrate through baryonic 
shells. The clump shells are assigned some angular 
momentum (similarly to the baryon and dark matter 
shells) that stops them from reaching the singularity 



spectivcly. Q l is the rate at which clump i heats the gas 
parcel it is embedded in at that time. We assume that 
many clumps are present within each parcel of gas, and 
that their motions are isotropic so their forces cancel 
out: 



Yf 1 



0. 



and 



E u 



1 "cl 



E^ 



(5) 



£tf (6) 



0. 



In the reference frame of a gas parcel , the clump 
loses energy, which is converted to heat, so using eq. (fTJ) 
we identify 



Q % = /dra 



ci = - c d a phot K 



(7) 



as the heating rate that clump i heats the gas par- 
cel in which it is embedded. The total energy of the 
clumps and of the gas should be conserved on aver- 
age according to eq. ([5]), assuming there are enough 
clumps so the averaging is correct, and that the assump- 
tion that the clumps have isotropic velocities is good. 
The difference equations in Hydra conserve energy algc- 
braicallj0(making energy conservation independent of 
resolution), and we do not want to violate this property 
We thus require a detailed balance between the clump 
deceleration and energy deposition of each clump: 



/>cl 



(8) 



2 an algebraic scheme is such that the exact energy term is always 
added to one component and subtracted from the other explicitly, 
making sure the energy is conserved to the machine accuracy 
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which implies: 

f — f Vrel 

Jc\ — /drag 

Vcl 

1 n A 2 l w rcl| ,„n 

= ~o CdAp hot v Tel (9) 

2 I'd 



Eq. recovers the exact solution when the gas 
parcel is at rest (which is the typical case of heating of 
a hydrostatic gas) and when there is no relative velocity 
between the gas and the clump (i.e. no drag). 

Physically, the drag forces always act to decrease 
the radial and tangential components of the velocity 
of clumps. While the radial velocity is replenished by 
the gravitational force, the tangential velocity monoton- 
ically decreases in time, so clumps lose angular momen- 
tum, becoming more radial in their trajectories. Since 
the angular momentum of the clump shells in our simu- 
lations is conserved, the clumps in the simulations can- 
not spiral towards the centre. To compensate for this 
problem, and allow clumps to reach the central core, 
a much smaller angular momentum is assigned to the 
clumps, placing them on almost radial trajectories. 



2.1.3 Fragmentation of Clumps 

Once the framework for accelerations and heating is de- 
fined, we proceed to implement recipes for clump evolu- 
tion. In this work, we implemented only the most cru- 
cial additional recipes from DB08: clump fragmenta- 
tion, and clump destruction. A clump fragments into 
^frag clumps (2 thorough out this work) once it has re- 
pelled its own mass of ambient gas. The amount of gas 
repelled is calculated by numerically integrating over 
T^ci /Ogas^rei dt in the same Runge-Kutta scheme used 
for the dark matter. After each timestep, the column 
mass and clump mass is compared. When the column 
mass exceeds m c \, m c \ is divided by rif rag , n c \ is mul- 
tiplied by nf ra g, and the column mass integral is reset 
to 0. As m c i becomes exceedingly small, its drag de- 
celeration become large, and the periods between con- 
sequent fragmentation events become shorter. Clumps 
are destroyed when their mass decreases below a critical 
mass, at which conduction and evaporation is expected 
to disintegrate them com pletely (in this work - 10 4 M Q , 
iDekel fc Birnboim|[2008l) . 



2.1.4 Destruction of Clumps 

The destruction of clump is achieved by adding its mass, 
to the corresponding baryonic shell. The velocity and 
angular momentum are not changed, and the internal 
energy and temperature is calculated by mixing the cold 
and hot components according to energy conservation 
by solving for the final internal energy, E? nt in: 

M (E int + E kin + Egrzv) + m (e int + e k i n + e grav ) = 
(M + m)(E kin + £ grav ) + (Af + to) E f mt . (10) 



with M , E correspond to the baryonic shell values, m, e 
to the clump shell. All energies are specific energies 
(per unit mass) and e^t = C v T c \ with T c \ = 10 4 K. 
Rarely, the final temperature will drop below 10 4 _ftT, at 
which case the final temperature is set to 10 4 K. This 
is found to occur only for baryonic shells that are cold 
(around 2 x lO 4 ^) and on a free fall to the BCG, and the 
temperature floor that is artificially applied never stops 
their infall. In this case, the fictitious energy is tracked 
throughout the run and is always negligible. The overall 
accreted mass is not effected by this correction. 

2.1.5 Applicability to 3D Simulations 

In the ID case described above, the gas parcel is a finite- 
dimension gas shell, and the clump is a thin shell. While 
3D simulations are beyond the scope of this work, if one 
wishes to incorporate the effects of cold clumps in 3D 
simulations, a generalization for the 3D case is readily 
available, as follows. Both SPH and grid-based (Eule- 
rian or Lagrangian) simulations that model dark matter 
as N-body particles can assign clumps to a dark mat- 
ter particle according to cq. ([2]), calculate its accelera- 
tion according to eq. © and heat the gas according to 
cq. ([7]). In cosmological simulations, it is also necessary 
to self-consistently create those clumps. These can ei- 
ther be crea ted semi-analytically according to cooling- 
instabilities (|Fieldlll965HBmnev et aTll2009ft or by iden- 
tifying unresolved ga seous clouds and replaci ng them 
with clump particles ([Keres k. Hernquistll2009L identify 
clump formation on Milky Way scales, but resolution 
would not allow to scale this procedure to galaxy clus- 
ter scales). 

2.2 Cell splitting - a ID Adaptive Mesh 
Refinement 

In the Lagrangian formalism, the amount of baryonic 
mass within each shell is constant throughout a sim- 
ulation. Clump destruction, however, violates that by 
dumping mass into the baryonic shells at the event of 
clump destruction. This mass deposition is expected to 
occur preferentially at specific regions. Indeed, the sim- 
ulations discussed below initially caused the formation 
of a few extremely massive shells that absorbed most 
of the mass of the clumps. This situation (of large con- 
trasts between adjacent shells in mass and widths) re- 
duces the accuracy of the numerical scheme, and ef- 
fects the accretion rates onto the BCG. To overcome 
this, an adaptive splitting of baryonic shells is per- 
formed: when the shell is wider by some factor than 
both its neighbours, or when the shell's width divided 
by its radius exceeds some preset fraction, the shell is 
split to two. The factors ultimately used were: Ar„ > 
4 x max(Ar n _i, Ar n+ i) and Ar n /r n > 0.2 (with Ar„ 
and r n the width, and central radius of shell n) which 
caused splitting to occur a few dozens times throughout 
a full, Hubble time simulation. This choice of parame- 
ters is a result of trial and error motivated by the re- 
quirement that except for transient periods, the width 
of Lagrangian shells are roughly equally spaced (lo- 
cally), ensuring the robustness of the Numerical scheme. 
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A shell is split into two constant mass shells. Values 
that are naturally defined at centres of shells (density, 
temperature) are treated as step functions and their 
values are equal in the two new shells. Values that are 
defined on boundaries (radius, velocity, angular momen- 
tum) are interpolated linearly with mass. This defini- 
tion ensures that the total internal energy of the system 
remains the same. The potential and kinetic energies, 
however, can change. This energy non-conservation is 
not corrected for explicitly. Rather, it is treated as non- 
conservation and tracked throughout the run. In the 
high resolution simulations shown below (2,000 bary- 
onic shells and 10,000 clump and dark matter shells), 
the total, overall energy is conserved to better than 10~ 2 
over a Hubble time. 



2.3 Convection and Mixing Length Theories 

2.3.1 Convection 

A long term balance between cooling and heating re- 
quires, in addition to sufficiently large energy injection, 
a correct distribution of this energy. The cooling rate, 
at constant pressure, scales as e coo i ~ p 1,5 , and most 
heating mechanisms (the drag forces discussed here, dy- 
namical friction, radiative heating from supcrnovae and 
AGNs) generally heat according to eh ca t ~ P- Even 
if at some point in space, and at some initial time, 
e C ooi/ehcat = 1, this ratio scales like p 05 . This rela- 
tion indicates a positive feedback so if density increases, 
cooling is more efficient, causing a further increase in 
density at constant pressure, and unstable cooling will 
occur. Vice versa, a density decrease increases the rela- 
tive importance of heating over cooling, decreasing the 
density further, and an over- heating instability occurs . 
This point has been made b y Conrov fc Ostrikeil (|2008ft 
(see, however, iKunz et alj|2010l who claim that resid- 
ual magnetic turbulence might inject energy into the 
gas in a stable manner) , using ID hydrodynamic cal- 
culations of clusters initially hydrostatic within a static 
potential well. They find that stable, long term equilib- 
rium requires fine tuning of the heating efficiency which 
is unlikely. The present work differs in that it treats the 
gravitational drag feedback and the cluster evolution 
from initial cosmological perturbation consistently, but 
should still suffer from a heating instability. Such an 
overheating will manifest as a shell or region of shells 
of gas continuously becoming hotter and under-dense, 
with very high entropy. In a 3D configuration, this en- 
tropy inversion is unstable to convection when entropy 
is declining outwards: a slightly under-dense parcel of 
gas floats buoyantly, carrying energy and momentum 
outwards, and over-dense parts sink towards the cen- 
tre reducing the average entropy and specific energy of 
the core. For our ID simulation we invoke a ID subgrid 
mode l for convection - mixing length theory ([Spiegell 
11963ft . 



2.3.2 Mixing Length Theory 

Convection occurs when an adiabatic displacement of 
a parcel of gas, in pressure equilibrium with its new 



position, results in a net force on that parcel tending 
to increase the displacement. This requires that, for an 
upward displacement, the temperature of the gas parcel 
will be smaller than its surrounding: 



AVT = 



dT\ _ fdT\ 
dr J \dr J s 
dT\ dS (dT 



dS 
dT 
dS 



P dr 
dS_ 

dr ' 



dP 



dP 
dr 



dT 
dP 



dP 
dr 

(11) 



with the first term in the r.h.s corresponding to the ac- 
tual temperature derivative in the profile, and the sec- 
ond to the adiabatic change in temperature as a result 
of the pressure profile of the halo. A negative value of 
AVT allows for convection to occur. It is convenient to 
derive this relation using a form of the ideal equation of 
state in which the two thermodynamic free parameters 
are the entropy, S, and the pressure, P. 

p N A k B „ 

P = pT, 



S 



P 

N A k E 



In 



T 3/2 



(12) 



/i p 

can be inverted (setting /t = p/N A &b) into: 

p=(AP) 3/5 e-( 2 / 5 ^ s , 

T= (/}P) 2/ V 2 / 5 ^ 5 . (13) 

Plugging these relation into the derivatives in eq. (JTTJ 
we get: 



2 as 

AVT= 5 MT ^ 



(14) 



recovering the known results that when the composi- 
tion of the gas is constant, entropy inversion leads to 
convection. 

The buoyant bubbles rise for a typical length be- 
fore being destroyed by Kelvin Helmholtz instability 
or conduction. The details of this destruction depend 
on the size of the bubbles, the smoothness of the den- 
sity and gravitational profile, conduction and magnetic 
fields. Solving for it requires fine 3D simulation of the 
convective process. We replace this dependency by a 
free dimensionlcss parameter, the mixing length (£), 
assuming that bubbles rise a distance, Hp , which is pro- 
portional to the atmospheric scale length of the halo at 
each point: 

H P =L—, (15) 

pg 

with g = GM/r , A typical acceleration (assuming iso- 
baric perturbations) is: 



Sp 



9— = 



9 

P P 

so a typical bubble velocity is 



dp 
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dr 



H P 
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(16) 



(17) 
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Once the velocity exceeds the local sound of speed in 
the halo, shocks are created which quickly act to mix 
the bubble with its surrounding. In the following calcu- 
lation the velocity is not allowed to exceed the speed of 
sound, c s , and in that case Hp is reduced until v = c s 
in eq. (fT7|) . A maximal convection model is a model 
with arbitrary high L, so effectively the bubbles always 
accelerate until the speed of sound, at which time they 
are broken and mixed. 

The Flux of energy per unit surface per unit mass 

is: 



F c = C P v [AVT Hp] =Cpv 



2 -T dS TT 
5 or 



(18) 



which is determined by the halo profile from the simu- 
lation at each time, and the mixing length L. Cp is the 
constant pressure heat capacity and is related to fx by 
C P = 5/(2£)._ 

A numerical solution of the mixing length model 
requires evaluation of the incoming and outgoing fluxes 
from the boundaries of each radial shell. The fluxes de- 
pend on the temperature gradient between each shell 
and the ones directly below and above it, and interpo- 
lation of thermodynamic properties from the shell cen- 
tres to the shell's edges is required. A solution using an 
explicit numerical scheme (with the fluxes determined 
at the beginning of each timestep) requires extremely 
small timesteps to avoid negative temperatures, so an 
implicit scheme which solves simultaneously for all the 
temperatures and fluxes at the end of each timestep in 
each convective area was implemented. This is done by 
inversion of the three-diagonal matrix which is obtained 
by discretization of eq. (|18p. 

In the context of the clump heating discussed in 
this paper, it should be emphasized that the convection 
is not between the cold gas in the clumps and the hot 
surrounding. It is only within the hot component, and is 
the result of heating and cooling of that component. The 
convection model assumes linear perturbations within a 
single gas phase which separates into two phases (hot 
buoyant bubbles, and cold sinking gas). The propaga- 
tion of clumps, which have a typical over-densities of 
XL 1000 is followed explicitly using the processes de- 
scribed in 32X21 - 32X41 

The ICM is mildly magnetized, with the non- 
thermal magnetic pre ssure contributing at most 10% 
of the total pressure (jChurazov et al.l 2008). This ef- 
fect could lead to heat-flux buoyant instability (HBI; 
iParrish et all I2009D even when the entropy profile is 
monotonically increasing provided there is a temper- 
ature inversion near the core (dT/dr > 0). These in- 
stabilities act to align the magnetic fields perpendicu- 
lar to the temperature gradi ent, in a manner th at sup- 
presses further conduction (jParrish et al.l [2009^ . Con- 
vective motions are also somewhat suppressed even for 
material that is hydrodynamically convectively unsta- 
ble (ds/dr < 0). Future work could, and should, use 
a revised mixing length theory for which the driver of 
convection is temperature inversion rather than entropy 
inversion. Such a model would need to take into account 
the saturation of the instability as the magnetic fields 
align themselves, baring in mind that heating by clumps 



is closely related to turbulence driving by clumps. Also, 
the l argely reduced convection strength that is eluded 
to in (jParrish et al.ll2008l [2009IPI must be evaluated. Ul- 
timately, the convection here is invoked to smooth over 
local instabilities for which, at least according to the 
spherical calculations, steep gradients of more than an 
order of magnitude in temperature and entropy form 
at the spatial resolution limit (thin red line of Fig. [3]). 
These extreme gradients (that are also present at edges 
of radio bubbles in clusters) are far from linear pertur- 
bations, and the validity of the linear analysis of the 
various convection prescriptions is highly questionable. 



2.4 The Observed Properties of Hot Gas with 
Cold Clumps 

In the current implementation, when clumps are de- 
stroyed, their mass is added to the hot component in- 
stantly. In reality, the process of KH fragmentation, 
followed by small scale evaporation and conduction of 
the debris will yield a multiphased gas, with an ef- 
fective entropy and temperature which is between the 
values of the hot and cold phases. As will be shown 
later (Q , the clump density and clump destruction 
rate increase towards the centre so effective cooler and 
lower entropy values are expected there. The radiative 
signature of gas heating through all the temperatures 
between 10 4 A" to the cluster ambient gas temperature 
of ~ 3 x 10 7 A is expected to be significantly differ- 
ent from that of radiative cooling since it is governed 
by heating processes (emission spectrum from heating 
gas, a lbeit by other heating mechanism s have been stud- 
ies in lVoit fc Donahudll997tlOhll2004D . A framework of 
heating and cooling processes in la yers between hot and 
cold media have been pro posed by iBegelman fc Fabianl 
(119901) : iGnat et al.l (12010ft . 

The observational signature of clump break up 
would require detailed 3D simulation of clump interac- 
tions with cluster core gas, and multiphased modeling 
of the radiative signature during the heating process, 
and is beyond the scope of this work. Instead, we will 
plot below the mass weighted entropy and temperature 
of the two components. This is a lower limit for the 
observed entropy and temperature as the clumps' con- 
tribution to the luminosity, particularly at X-ray wave- 
length, is probably small. Physically, it corresponds to 
the thermodynamic properties expected in the event of 
full mixing between the cold and hot phase. The actual 
temperature and profile expected from the multiphase 
gas is thus bracketed between the hot only component, 
and the mass weighting between the hot and cold com- 
ponents presented in figs. [3] and |3J 



3 RESULTS 

The simplified model described in fj2] spans a multidi- 
mensional parameter space including the fraction of ac- 



3 See, for example , the short-dashed lines in fig. 6 and 10 of 
IParrish et all fcOQgft , 
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creted gas in clumps, the initial clump masses, the num- 
ber of fragments that a clump breaks up to, and the 
mixing-length parameter for convection. In the absence 
of additional physical insight concerning the formation 
mechanism and properties of these clumps, we attempt 
to find a working set of values for the model parameters, 
to serve as a feasibility test and hopefully provide clues 
for acceptable values of the key parameters. Ultimately, 
a more systematic survey of parameter space will have 
to be conducted, with physically motivated values for 
key parameters such as number and mass of clumps. 

In this section we restrict ourselves to one typi- 
cal CC cluster halo with virial mass of 3 x 10 14 M Q by 
z = 0, a diffuse baryon fraction of 10%, and a smooth 
accretion history according to the average g rowth rate 
of the main progenit or a la iNeistein et al.l (|2006T > (see 
iBirnboim et al.l [20071 for a detailed description). The 
metallicity is assumed to be constant at Z = 0.3 Zq, 
and the c ooling is Bremsstrahlung and line cooling ac- 
cording to ISutherland fc Dopital (|1993|) . The initial res- 
olution is 2000 baryonic shells and 10,000 dark-matter 
shells, roughly logarithmically spaced in their initial 
radii. When shells expand, the adaptive mesh refine- 
ment algorithm splits them ( N2.2p , so the resolution near 
the cluster core (~ 50kpc) at all times is better than 
~ 2kpc. This yields converged results in terms of the 
profiles and the amount of gas that cools. We implement 
three different models for clump heating, as follows: 

(i) Model C is the null model with no clump heating 
and no convection. It is meant to to reproduce the over- 
cooling problem. 

(ii) Model CH adds clump heating but no convection, 
so it should show the over-heating instability. The frac- 
tion of baryons in clumps is 5% and the clump initial 
mass is 1O 8 M . Th e clumps are simulated by 10, 000 
clump shells ( i|2.1.2|) . 

(iii) Model CHC has the same clump heating as in 
CH but with maximum convection turned on f i|2.3.2[) . 
The smoothing by convection is supposed to regulate 
the clump heating and yield relaxed clusters compatible 
with observations. 

Figure [1] shows the time evolution of the gas in our 
simulated cluster comparing the three different mod- 
els for cooling, clump heating and convection. The ini- 
tial Hubble expansion and consequent turnaround of 
the Lagrangian gas shells is clearly seen, and the virial 
shock can be easily identified after a collapse by a fac- 
tor of ~ 2, both by a jump in temperature from below 
10 4 K to above 10 7 K, and by the abrupt slow down of 
the infall velocity, which is almost brought to a halt 
behind the shock. The global large-scale properties of 
the cluster are not affected by the addition of heating 
and convection. The virial radius evolves in a similar 
way, and the typical temperature in the halo at z = 
remains at T - 2 - 3 x 10 7 K (-2-3 keV), consistent 
with the expected virial temperature of 2.2 x 10 7 K for a 
cluster of virial mass 3 x 10 14 M Q . However, the models 
differ at the core, within the innermost 100 kpc espe- 
cially during the last 6 Gyr of evolution. Model C shows 
inwar d cooling flows at all times, as expected (jFabianl 
Il994t ). With the addition of clump heating in Model CH 
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Figure 2. Time evolution of X-ray luminosity (top), BCG mass 
(middle), and BCG accretion rate (bottom), all smoothed over 
lGyr. 



The cooling flows are stopped before t = —6 Gyr, the 
gas in certain shells is over-heated to extreme temper- 
atures £ 10 8 K, these shells are interlaced with cooler 
shells of — 10 6 K, and together the whole core inflates. 
This behavior is in conflict with the relaxed nature and 
smooth entropy and te mperature profiles of CC clusters 
(|Donahue et al.l r2006f ). The addition of convection in 
Model CHC removes the local over-heating, and keeps 
the core in equilibrium at the virial temperature with no 
cooling flow. A more detailed comparison of the models 
and observations follows. 



3.1 Time Evolution of X-ray Luminosities and 
BCG Masses 

Figure [5] shows the evolution of X-ray luminosity, the 
mass of the BCG, and the accretion rate onto it, for 
the three different models. The BCG is represented by 
the mass in the "disc", the mass that is supported by 
angular momentum in the ID simulations, extending to 
— 10 kpc. The X-ray luminosity is obtained as the total 
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Figure 1. The time evolution of the gas in a cluster for the three different models of clump heating. Shown in grey lines are the evolving 
radii of 25 spheres that encompass fixed baryon masses equally log-spaced in the range 10 10 — 5 X 10 13 Mq. The dark matter and clumps 
are not shown. The final halo mass is 3 X f0 14 MQ. Look-back time and redshift are marked at the bottom and at the top. The color 
refers to gas temperature. Top: Model C, cooling only. Middle: Model CH, cooling and heating. Bottom: Model CHC, cooling, heating 
and convection. Model C shows cooling flows in the core between 10 and 100 kpc during the last 6 Gyr (in particular near t ~ —2 Gyr). 
Model CH shows core over-heating and expansion in the last 6 Gyr, and Model CHC demonstrates how convection regulates the heating 
and brings the cluster to an equilibrium. 



cooling radiation from the gas outside the BCG0. The 
quantities are smoothed over 1 Gyr to erase sharp fea- 
tures that result from the discreteness of the calculation 
in the idealized spherical calculation. 

Model C shows a variable luminosity, which occa- 
sionally exceeds by an order of magnitude or more the 
observed luminosity of 1 44 — 10 45 erg sec" 1 as derived 
from the L x — T relatio n (jEdge et alJll990l iDavid et all 
[19931 lMarkevitcljL99"8l) . In Model C the final BCG mass 
exceeds 2 x 10 12 M Q and the accretion rate, which is an 
indicator for the star formation rate, has long episodes 
where it is in the range 100—1000 Mq yr _1 , at odds with 

4 Most of the energy is emitted from cooling of ~ keV gas, and 
is predicted to contribute to the X-ra y luminosity that is relevant 
to the L x — T relation (for example: [Markevitch][l99^). 



observed CC clusters. When clump heating is added in 
Model CH, the cooling flow into the BCG is drasti- 
cally suppressed since before t ~ — 8 Gyr, and it reaches 
a complete shutdown after t ~ — 5 Gyr. The luminos- 
ity maintains a high level of ~ 2 x 10 45 erg sec -1 since 
t r-j —7 Gyr. This results from the over-heating insta- 
bility in the halo core, leading to very dense shells that 
boost the dissipation due to drag interaction with the 
clumps and enhance the resulting radiation. The ad- 
dition of convection in Model CHC brings the BCG 
mass to a constant value of ~ 10 Mq with no de- 
tectable cooling flow since t ~ —6 Gyr. The smoothing 
of the instability brings the luminosity to a low value of 
~ 10 44 erg sec -1 , consistent with the observed Lx — T 
relation. 
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Figure 3. Radial profiles of entropy at z = for models C, 
CH and CHC. Shown in three curves is the entropy of the hot 
medium, and an additional curve marked CHC e gf refers to the 
effective entropy of the mixture of hot gas and the surviving cold 
clumps, computed as a mass- weighted average and brackets the 
predicted entropy from below. 



3.2 Entropy and Temperature Profiles 

The z = entropy profiles of the three models are plot- 
ted in Fig. [3] and the temperature profile of model CHC 
is plotted in Fig. @] The entropy profiles in the outer 
halo, overall increasing close to linearly with radius and 
having values ~ 100 keV x cm 2 at 100 kpc, are similar 
in all m odels and consistent with observations of CC 
clusters ([Donahue et all 120061 : ICavagnolo etHI I2009T ). 
The apparent periodic fluctuations represent cold fronts 
that result from mergers of outward-propagating shocks 
at regular intervals with the virial sh ock (also visible 
in Fig . [T]), which have been studied in iBirnboim et all 
(|2010T> . While the fluctuations in Models C and CH 
are locally non-monotonic, the entropy profile of Model 
CHC is monotonically increasing throughout because 
the convection removes negative entropy gradients. 
Model CH shows a high-entropy core due to overheat- 
ing, and a strong variability representing a mixture of 
cold-dense and hot-dilute shells that result from the 
overheating instability, both in conflict with observa- 
tions. The convection introduced in Model CHC re- 
moves the fluctuations and produces what seems to be a 
flat core entropy profile inside 150 kpc at 50keV x cm 2 . 
Such a core is still inconsistent with CC cluster cores, 
but recall that the profile shown is limited to the en- 
tropy of the hot component alone. The effective entropy 
that could actually be observed is a mixture of the en- 
tropy in the hot gas, cold clumps, and anything in be- 
tween as discussed in §2.41 The effective entropy profile 
shown in Fig. [3] is monotonically increasing in the core 
down to ~ 30 kpc, consistent with CC clusters. 

The temperature profile of Model CHC at z = is 




log R [kpc] 

Figure 4. Radial profiles of temperature at z = for model 
CHC, for the hot medium and for the mixture of hot medium and 
cold clumps, computed as a mass-weighted average. The predicted 
temperature is bracketed between this upper and lower limit. 

plotted in Fig. @] It shows a roughly isothermal halo at 
the virial temperature outside the core of ~ 100 kpc, 
with a mild decline toward the virial radius, as ob- 
served ([Donahue et al.|[2006() . This large-scale tempera- 
ture profile has not been affected much by clump heat- 
ing and convection. The temperature of the hot compo- 
nent is rising toward the centre, by a factor of ~ 2, but 
the effective mass- weighted temperature of the mixture 
of hot and cold components is declining toward the cen- 
tre, by a factor of a few. This mass- weighted effective 
temperature (that corresponds to the single tempera- 
ture the gas would if the phases were fully mixed) is a 
lower limit to the luminosity- weighted temperatures ob- 
served, and is consistent with the moderate temperature 
decline in CC cluster cores. The temperature profiles of 
the various models, and its time evolution, can also be 
seen in Fig. [1] 



3.3 Sensitivity to choice of parameters 

The results presented in this section indicate that 
clumps of m c = 10 8 M Q that make for f c = 0.05 
of the gaseous component of clusters is sufficient to 
remedy the over-cooling problem, and quench the ac- 
cretion of gas onto the BCG. These parameters were 
picked to demonstrate the effectiveness of our model. 
We find that our simulated clusters are not particu- 
larly sensitive to these values, and that no fine tun- 
ing is required. Rather, a wide envelope of allowed pa- 
rameters is allowed. Fig. [5] is analogue to Fig. [5] and 
compares the luminosity, BCG mass and accretion rate 
predicted by simulations with maximal convection for 
a parameter choice of: (m c = 10 7 M Q ,/ C = 0.1) and 
(m c = 1O 8 M0, f c = 0.1), along with our fiducial model 
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Figure 5. Time evolution of X-ray luminosity (top), BCG mass 
(middle), and BCG accretion rate (bottom), all smoothed over 
1 Gyr of three sets of parameters. 



of (m c = 10 8 M Q , f c = 0.05). The same behaviour is also 
found in the resulting radial profiles. A simulation with 
the parameters (m c = 1O 7 M ,/ C = 0.05) allowed for 
too much cooling, and final BCG mass of 1.5 x 1O 12 M 
- that is slightly excessive. 

Beside these parameters, some of the theoretical 
model assumptions (for example the drag efficiency, 
the fragmentation of clumps and the convection) might 
need to be modified after more detailed 3D simulations 
or additional observational constraints are found. We 
expect that for a different model the parameters will 
need to be readjusted, but since the heating model is 
not particularly sensitive, we expect that such a choice 
will always be possible. We do not believe that a com- 
prehensive parameter survey will be beneficial at this 
point, until the different components of this model are 
better constraint either theoretically or observationally. 




1 2 3 

R [kpc] 

Figure 6. Average number density of clumps (top), average 
clump mass (middle) and fractional volume in clumps (bottom) 
as a function of radius at z = for Models CH and CHC. 



4 OBSERVATIONAL SIGNATURES OF 
CLUMP HEATING 

The model parameters chosen in this paper were mo- 
tivated by the analysis of DB08 and were calibrated 
such that the model crudely reproduces the properties 
of CC clusters in order to demonstrate the feasibility 
of such a model of gravitational heating. The repro- 
duced properties include the BCG mass, the cold mass 
accretion history, the X-ray luminosity and the entropy 
and temperature profiles. This model makes additional 
predictions that could distinguish it from other heating 
models such as the ones based on AGN feedback. Some 
of these predictions are discussed here. 



4.1 Cold gas in the ICM 

The mass function of clumps at every radius is a dis- 
tinctive prediction of our proposed model. Assuming an 
initially uniform population of 10 8 A/ Q clumps, Fig. [5] 
shows the profiles of number density and average clump 
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mass for Models CH and CHC. The fraction of volume 
occupied by cold clumps is shown in the bottom panel. 
These predictions depend on the specific choice of the 
initial mass of the clumps and on the baryonic fraction 
of mass in these clumps. The number density of clumps 



increases toward the centre oc 



partly reflecting 



the general density profile of the cluster and partly be- 
cause of the break up of clumps into fragments f ^2.ip . 
Once clumps are fully destroyed (fragments < 10 4 M Q ), 
their mass is added to the hot component, and they 
are no longer plotted in Fig. [6j The figure thus shows 
the steady state population of clumps as they are con- 
tinuously accreted and destroyed. The average clump 
mass is declining from the initial value of 1O 8 M in 
the outer halo to ~ 10 6 M Q near the centre, reflecting 
the clump fragmentation as they flow in. The fraction 
of volume occupied by the cold clumps peaks at a few 
percent near the cluster core, and drops to smaller val- 
ues at larger radii. This justifies ignoring the additional 
pressure caused by this component. 

The clumps are initially cold, and, except for the 
mild compression they undergo as they fall in following 
the increasing pressure of the ambient gas, they are not 
expected to heat up or emit much radiation. However, as 
the clumps are disrupted by hydrodynamical instabili- 
ties and possibly also by tidal effects, they fragment into 
smaller pieces for which conduction and evaporation be- 
come more important ( ij2.4[) . Once heated to interme- 
diate temperatures, the gas begins to radiate. Spectro- 
scopic observations could in principle constrain the va- 
lidity of this model in comparison with AGN-feedback 
models, where one expects gas cooling rather than heat- 
ing through the intermediate temperatures. The cur- 
rent simplified implementation of clump heating does 
not permit a proper comparison, which is left for fu- 
ture work. Additionally, three-dimensional simulations 
are required for a detailed analysis of the shape of the 
clumps as they are stretched pe rhaps leading to mor - 
phologies resembling filaments (|Murrav fc Linl I2004T ) . 
This emission, in H a and line and continuum emis- 
sion of the intermittent X-ray temperature gas, may 
allow more accurate comparisons of this model with 
the observed profiles in cluster c ores. Observations in 
the Perseus cluste r (NGC 1275) (|Conselice et alj|200lt 
iFabian et al.ll2008T) show a complicated structure of H a 
filaments and blobs. The typical masses of these fea- 
tures are 10 6 — 10 8 M Q , consistent with the allowed mass 
range for clumps in DB08, and with the distribution 
predicted by our model (Fig. |6]). We note that this re- 
sult depends on the initial mass of the clumps - a free 
parameter here. The consistency of this prediction with 
observation is an indication t hat our choic e of in itial 
mass of lO 8 A/0 is reasonable. IFabian et al.l (|2008l ) in- 
voked strong magnetic fields to stabilize the filaments 
for cosmological times, such that their age can match 
that of the observed radio bubbles. Our heating model 
suggests instead that these filaments are constantly be- 
ing destroyed, as new clumps enter the cluster core, get 
stretched and destroyed, and create new filaments. The 
projected filling factor of these structures approaches 
unity within the innermost 10 kpc and it drops out- 
wards (jConselice et al.|[200l[ ). Such a behaviour is pre- 
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Figure 7. The average distance between clumps (top), the aver- 
age drag energy injection compared to the cooling rate (middle) 
and the ratio of turbulent pressure to thermal pressure (bottom) 
all as a function of radius at z = for Model CHC. 



dieted by our model (Fig. [S]). H n emission in othe r clus- 
ters have been reported bv lHeckman et all (|1989D . who 
found that the cold gas has velocities at random di- 
rections rather than a coherent radial cooling-flow pat- 
tern. This kinematics could be interpreted as clumps 
oscillating in and out at the vicinity of the BCG. Struc- 
tures of neutral gas are also see n in the Virgo Cluster 
by the ALFALFA 21cm survey (|Giovanelli et al.l 120071; 
iKent et aLI |2007( ) showing evidence for neutral gas ar- 
ranged in clumps, with masses as low as the detection 
limit of 2 x 10 7 Mq, sometimes with no optical counter- 
parts. 



4.2 Turbulence in the ICM 

Another prediction of this model is the power of tur- 
bulence that is produced in the ICM. When the clump 
velocities are subsonic with respect to the ICM, the en- 
ergy and momentum of the drag are first converted to 
kinetic energy in turbulence, which cascades down to 
smaller scalelengths where it dissipates into heat. When 
the motion is supersonic, some of the energy is con- 
verted directly into heat, but since momentum is con- 
served, some of the energy must be transferred as kinetic 
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energy to the ICM. We assume that the turbulence is 
generated on a scalelength comparable to the average 
distance between clumps, L, and that a Kolmogorov 
spectrum governs the cascade of eddies from this scale 
to smaller scales. The turbulent energy and pressure are 
then 



eturb — ^ C o( e drag L) 2 ^ 3 , 



(19) 



-Pturb = (7 — 1) PICM eturb = C PlCM^drag L) 2 ^ 3 

(lLandau fc Lifshit3 Il959l ). with cq = 2.1 taken from 



iPopol ([20001 )7" The amplitude of the turbulence spec- 
trum, and the turbulent pressure can thus be deter- 
mined once we know the average distance between 
clumps and the drag heating rate at every radius. The 
ratio between the turbulent pressure and the thermal 
pressure for Model CHC is plotted in Fig. [7] The level 
of turbulent pressure in the cluster core is at the level of 
~ 6% of the therm al pressure, in qualitative agreement 
with observations (jRebusco et al.|[200~ ; IChurazov et al.l 
12008ft . The average clump distance, L and the energy in- 
jection rate ed rag are separately shown is Fig. [7] as well. 
The value of L at the center is i n rough agreeme n t with 
the H a structures observed in IConselice et al.l (|200lf ) 
(further discussed below) and the drag heating exceeds 
the cooling rate near the centre, but becomes compara- 
ble around the 100 kpc, consistent with convection oc- 
curring near the centre as suggested in 



4.3 High Velocity Clouds in the Galactic Halo 

High velocity clouds (H VCs) are observed in 21cm HI 
data ([Blitz et al.lll999t ) as concentrations of gas mov- 
ing at velocities > 100km s -1 relative to the rotating 
frame of the Milky Way. Options for the spatial ori- 
gin of the HVCs ran ge from the Milky Way (MW) disc 
(jWakker et al.l [20081 an d references therein), through 
the M agellanic clouds (lOlanol 120081 ). to extragalactic 
origin (|Blitz et al.lll999D . The distance to and the ion- 
ization fraction of individual c louds (and therefore t heir 
size and mass) are unknown. iPutman et al.l (|2003[) es- 
timated distances of HVCs based on their H Q flux and 
models for the emission of ionizing radiation from the 
MW. They find, within the modeling and measurement 
uncertainties, that most HVCs are within a distance 
of < 30 kpc, indicating a mass range of 10 4 — 10 8 M Q , 



n dicatmg 

(|Putman et alj 120031: iBirnboim fe Loebl I2009D Other 



estimates by absorption features ([Thorn et al.l l2008f) 
yield comparable distances and masses. The origin of 
the HVCs is unclear. Models have proposed that they 
form within the Galactic halo by cooling instabil i- 
ties dMaller fc Bullockl 12004 IKeres fc Hernquistll2009D 



See, h owever jFraternali fc BinnevT i 20081 ) . Binnev et al.l 
(|2009|) analyzed formation of clumps in smooth MW 
halo conditions from thermal instability and concluded 
that it is unlikely to occur near the centre. They find, 
however, that the conditions become favorable closer 
to the halo virial radius and when th e entropy profile 
is sha llow (as is shown numerically in iKaufmann et all 
l2009f) . Their stability analysis tests growth of instability 
from infinitesimal perturbations but the growth of non- 
linear perturbations caused by shocks, collisions, and 



gravitational perturbers depends on initial conditions. 
The line-emission peak of the cooling curve at ~ 10 5 A' 
would make the warm cosmic filaments outside clusters 
more susceptible for clump formation. However, clumps 
have been shown to be a natural consequence of cold- 
flow filament breakup by hydrodynamic instabilities 
([Keres fc Hernquist] [2009h . For example, streams that 
do not flow radially to the halo centre are susceptible to 
Raylcigh- Taylor instability. Also, shocks that originate 
from the galaxy, e.g., by mergers or by starbusts, arc 
likely to form clumps by Richtmycr-Mcshkov instability. 
The destruction of these clumps, and their interaction 
with the Galaxy and the IGM, are likewise under deba te 
(jFraternali fc Binnevl [20081: IKeres fc Hernquist]l2009D . 

We point at an obvious analogy between the ob- 
served Galactic HVCs and the clumps addressed in this 
paper. Their masses are comparable, and their spatial 
distributions in the halo and toward its centre are pos- 
sibly similar. The larger pressure in the ICM of a more 
massive halo would make the cluster clumps denser than 
the Galactic HVCs (by the ratio of virial temperatures 
which is more than 10 times larger for clusters), so 
clumps of a similar mass could survive longer in clus- 
ters, allowing them to travel to the centre according 
to the estimates in §5] The total mass encompassed in 
HVCs seems to be > 10 9 M Q , making it a few percent 
of the total baryons in the MW halo, in good agree- 
ment with our fiducial choice of parameters. A missing 
piece of the model is the yet unspecified origin of the 
clumps, in terms of physical mechanism and location. 
The existence of HVCs provides circumstantial evidence 
that such clumps might form. As long as the clumps 
are formed before they fall into the halo, or even if 
they form inside the halo at a radius that is not much 
smaller than the virial r adius (|Maller fc Bullockl 12004 
IKeres fc Hernquistll2009| ) , the gravitational energy that 
is released during their infall is significantly larger than 
the energy required for heating the clumps to the virial 
temperature (DB08). 



5 DISCUSSION AND CONCLUSION 

The concept of gravitational heating of ICM g as as a 
partial or full solution to the cooling-flow problem is 
more general than the specific clump model discussed in 
this paper. It is easy to show that the gravitational en- 
ergy that is released as baryonic matter falls in through 
the halo potential well is enough to balance the cooling 
rates in groups within haloes of virial masses ~ 10 13 M Q , 
and it exceeds the cooling rate by more than an order 
of mag nitude in cluster haloes ~ 10 14 - 15 Mm ( DB08) . 
This point h as also b een made in Fabianl (12005 



Wang fc Abell ([20081): iKhochfar fc Ostrikerl (| 20Qg 
El-Zant et aTl (12004): fFaltenbacher fc Mathewsl ([2007 
and Khochfar fc Ostrikerl ([20081 ) tap into the same en- 
ergy sour ce, but utili z e dyn a mical friction that is less 
effective. iNaab et al.l (|2007l) ; Uohansson et al.l ([2009D 
show, however, that dynamical friction can efficiently 
stop gas accretion onto massive elliptical galaxies. Con- 
duction also taps into this e nergy source, and was 
proposed three decades ago ([Bertschinger fc Meiksinl 
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119861: iRosner fe Tuckerl Il989h . It is probably ruled 
out beca use of magnetic f i eld suppression of th e con- 
duction (jBinnev fe Cowid Il98lt iFabianl 11994 and 
reference within), and because th e resulting profiles 
would have a flat te mperature core (iBregman fe Davidl 
19881). See, however. iNaravan fc Medvedevl ([20011 ) and 



Kim fc NaravarJ (feOOa T for alternative ideas. Here we 



use the same energy source, but the physical process 
used to couple the energy with the baryonic cooling 
component is hydrodynamic drag for which the strength 
of interaction peaks at low clump masses and high ve- 
locities. 

The challenge for every heating model is to dis- 
tribute the energy uniformly t hroughout the cluster 
core both in space and time (|De Young et al.l 120081 : 
ICattaneo et al.l [20091) . In order to obey the observa- 
tional constraints, the heating mechanism should sup- 
press the gas mass that actually cools by two orders of 
magnitude. Such a shutdown requires that the mecha- 
nism should act smoothly over a scale of a few kpc, set 
by the smallest object that would cool in the absence of 
feedback while being continuously heated by conduction 
from its surrounding. This scale follows from 



icond ~ a/ V «s P t = l\j rjo.2 r 2 5/2 n_l t 8 kpc 



(20) 



with Ksp the Spitzer coefficient (|Spitzerlll962f ) , rj the re- 
duction of the Spitzer coefficient due to magnetic fields, 
t the cooling time of the gas, and T 2 = kT/2 keV, n_ 2 = 
n/10~ 2 cm" 3 , t s = i/10 8 yr, 770.2 = rj/0.2. The value 
77 = 0.2 is a r easonable upper limit for t he efficiency 
of conduction (|Naravan fc Medvedevl l200lt ). It also has 
to be temporally smooth over the cooling timescale , 
which is a few 10 s yr at most dDonahue et al.l [20061 ). 
The fact that the energy source of gravitational infall is 
automatically distributed over the cluster volume and 
over Hubble times makes it easier to meet this challenge 
with gravitational heating than it is with AGN-feedback 
models, where the source is on scales ten orders of mag- 
nitude smaller than the cluster scales. However, the cou- 
pling of the infalling baryons with the ICM should be 
such that most of the gravitational energy is deposited 
in the cluster core. 

This paper addressed one specific scenario of grav- 
itational heating, in which the mechanism for feeding 
the energy into the ambient ICM is via hydrodynamic 
drag acting on small clumps of cold gas. In a typical 
cluster, with a modest gas fraction ~ 5% of the ac- 
creted baryons in cold clumps of ~ 1O 8 M0, the clump 
heating suppresses the cooling flows toward the BCG 
for the last 6 Gyr. The conditions at the core do not af- 
fect the incoming clumps, so the process is not strictly 
self-regulating. Regardless, we find that the large-scale 
properties of the cluster, such as the overall gas frac- 
tion, virial shock and temperature outside the core, 
are unaffected by the heating. Furthermore, the core 
does not explode; it reacts to the heating smoothly and 
quiescently without a need for inherent self regulation. 
Due to the over-heating, the central density decreases 
and the total X-ray luminosity emitted from the core 
(Fig. [2]) declines such that the core obeys the observed 
L x —T relation. On cluster core scales, convection acts 



to flatten the entropy profile of the hot component and 
carry heat outwards as expected. The effective entropy 
profile after taking into account the cold clumps as well 
does not exhibit an entropy core, and is consistent with 
observed entropy profiles, to within the model limita- 
tions that are discussed. A local instability caused by 
the linear dependence of the heating rate on density 
(eheat ~ p) acts to create extreme entropy and temper- 
ature peaks of sub-kpc scales. These peaks create strong 
convection that, once accounted for in the simulations, 
stabilize the heating process on local scales as well. 

With the fiducial values of the parameters used in 
this work in a 3 x 10 14 A/ Q cluster halo, the model is 
successful in quenching the cooling flows and in repro- 
ducing adequate BCGs, X-ray luminosities, and entropy 
and temperature profiles. The model also predicts the 
expected level of turbulence in clusters, and the fraction 
of cold g function of radius. These two observ- 

ablcs are predictions of this model, while they are not 
naturally addressed by AGN-feedback models. 

In DB08 we analyzed clump heating in a static halo 
using a Monte-Carlo approach to simulate an ensemble 
of clump trajectories. We realized that the heating rate 
in the core is higher than the cooling rate, which could 
cause the core to expand. In order to see how the clus- 
ter could reach a steady-state configuration one must 
allow the system to respond dynamically. The current 
implementation of the model using a ID hydro code 
to simulate a cluster in the cosmological context allows 
us to do just that. We find that convection can reg- 
ulate the over-heating instability and produce a clus- 
ter with no cooling flow in steady state. The net ef- 
fect of the dynamical response is to make the heating 
more efficient. For example, with 5% of the baryons 
in 10 8 Mq clumps inflowing into a static 3 x 10 14 M Q 
cluster, our estimates in DB08 indicates a heating to 
cooling ratio slightly below unity, while here we find it 
to be above unity, predominantly due to the net expan- 
sion of the core. The dynamical evolution of the cluster 
(Fig. Q]) is noticeable especially from z ~ 2 and on, as 
the core takes a different thermodynamic trajectory in 
response to the heating. The BCG mass is smaller, and 
the core density is lower than in the simulation with- 
out heating. It is therefore possible that other heating 
mechanisms that were tested within a static framework 
(IConrov fc Ostrikerl200ilFabiadl2003tlKim fc Naravanl 
120031: iKim et al.l 120051: ICiotti fc Ostrikerl |2007l ) would 
also show different evolution tracks once the gas is al- 
lowed to dynamically adjust to the energy input while 
the halo is growing. 

The main missing piece in the proposed model 
of clumps as the agents for depositing the gravita- 
tional energy of infall in the ICM core is the un- 
specified mechanism and birthplace for the formation 
of clumps with the desired properties of abundance 
and mass. While clumps p robably do not for m in- 
situ at the cores of clusters ([Binnev et all 120091) they 
can form within cosmic filaments (Dekel et al.l 120091: 
ICeverino et al.|[20Tot iFumagalli et al.ll201ll) and in the 
edges of haloes. Non-linear perturbations and compli- 
cated halo ge ometries can stimulate th e formation of 
such clumps ([Keres fc Hernquistl I2009D . The analogy 
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between the observed HVCs and the desired clumps for 
heating clusters is promising and may provide a clue for 
the origin of these clumps. 

The degree of dumpiness needed for effective clump 
heating, at the level of ~ 5%, does not seem to be very 
demanding. The clumps may be hard to detect out- 
side the halo virial radius as their temperatures are 
expected to be only slightly lower than the tempera- 
ture of the surrounding filaments, and therefore their 
inner densities in pressure equilibrium are expected 
to be only slightly higher. The clumps become denser 
and possibly more detectable once they enter the hot- 
ter and denser ICM, and especially as they approach 
the cluster core. We are encouraged by observations of 
H a structures with masses of 10 6 — 10 8 M Q around the 
Perseus BCG (|Fabian et all I2008D but a more careful 
comparison needs to be made, addressing the ioniza- 
tion states and the radiative signature of clumps as they 
are disrupted an d heated (jBegelman fe Fabianl Il990t 
iGnat et al.l l2010h . One should also work out the spec- 
trum of the X-ray emission from the multi-phased gas as 
it is heated by conduction and radiation. Here we only 
provide upper and lower limits to the observed temper- 
ature and entropy, derived by assuming either that only 
the hot component is observed, and that full mixing oc- 
curs instantaneously as the clumps disintegrate. 

The exact details of clump-ICM interactions, and 
the response of the ICM, cannot be properly addressed 
in ID simulations where hydrodynamic instabilities 
are almost completely suppressed independent of the 
quenching mechanism. In this first crude study we rely 
on the fact that the model proposed here naturally de- 
posits the energy over ~ 1 kpc scales (Fig. [5]) in a contin- 
uous manner. Nevertheless, a realistic study of whether 
the gas cooling is sufficiently suppressed would require a 
proper 3D simulation with clump heating implemented. 

The inherent runaway expansion that occurs when 
heating is faster than cooling is damped by a ID convec- 
tion model, which introduces a free mixing-length pa- 
rameter that can only be calibrated by 3D simulations. 
Our assumption here, that the convection is maximal in 
the sense that bubbles accelerate until they reach the 
speed of sound may be overly optimistic. Additionally, 
weak magnetic fields in the ICM may affect the na- 
ture and strength of the convection and could alter the 
general behaviour of convection to follow temperature 
inversions rather than entropy inversions. We find that 
the results are not particularly sensitive to the value of 
the mixing length parameter as the local perturbations 
are short scaled and (at least within the framework of 
the model) entropy inversion is erased for a wide range 
of mixing length parameters. These effects must be ad- 
dressed in future work. 

In A2.1.5I we outlined a proposed implementation 
of a 3D subgrid model for these clumps, which will al- 
low all these issues to be addressed. A different kind 
of 2D and 3D simulations, of interactions between sin- 
gle clum ps with the ICM gas , have been conducted in 
the past (|Murrav fc Linl |2004[ ) . but not for the specific 
conditions in clusters, and without some of the crucial 
physical components such as cooling and conduction. 

Our results support the notion that gravitational 



heating by the instreaming baryons could be a major 
player in the heating of the cores of massive galaxies and 
clusters. Whether this mechanism by itself is sufficient 
for preventing cooling flows or it must work in concert 
with AGN feedback is yet to be investigated using 3D 
cluster simulations (e.g. Zinger et al., in preparation). 
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